
from math import sqrt, asin, pi

def solve(r):
    # find all stations by for all <= i <= j <= k
    r2 = r * r
    for x in range(0, int(sqrt(r2 / 3.0)) + 1):
        x2 = x * x
        for y in range(x, int(sqrt((r2 - x2) / 2.0)) + 1):
            z2 = r2 - x2 - y * y
            z = int(sqrt(z2))
            if z * z == z2:
                yield x, y, z

class cached_risk(object):
    def __init__(self, r):
        self.__r = r
        self.__risk_dict = {}
    
    def __call__(self, a, b):
        dx, dy, dz = a[0] - b[0], a[1] - b[1], a[2] - b[2]
        d2 = dx * dx + dy * dy + dz * dz
        r = self.__risk_dict.get(d2)
        if r is None:
            r = self.__distance2risk(d2)
            self.__risk_dict[d2] = r
            return r
        else:
            return r
        
    def __distance2risk(self, d2):
        angle = asin(sqrt(d2) / 2 / self.__r) / pi * 2
        return angle * angle

class noncached_risk(object):
    def __init__(self, r):
        self.__r = r
    
    def __call__(self, a, b):
        dx, dy, dz = a[0] - b[0], a[1] - b[1], a[2] - b[2]
        d2 = dx * dx + dy * dy + dz * dz
        angle = asin(sqrt(d2) / 2 / self.__r) / pi * 2
        return angle * angle
    
def moon_dijkstra(r):
    # init stations
    stations = []
    for i, j, k in solve(r):
        ex = set() # use set to de-dupe
        ex.add((i, j, k))
        ex.add((i, k, j))
        ex.add((j, k, i))
        stations += list(ex)
    # print r, len(stations), stations
    # dijkstra
    risk_table = {(0, 0, r): 0}
    # risk_calc = cached_risk(r)
    risk_calc = noncached_risk(r)
    current = (0, 0, r)
    while len(stations) > 0:
        current_risk = risk_table[current]
        stations.remove(current)
        min_risk = r
        min_station = None
        for station in stations:
            risk = risk_table.get(station, r)
            new_risk = current_risk + risk_calc(current, station)
            if new_risk < risk:
                risk = new_risk
                risk_table[station] = new_risk
            if risk < min_risk:
                min_risk = risk
                min_station = station
        current = min_station
    min_journey = risk_table[(0, r, 0)] * 2
    for s in risk_table:
        risk = risk_table[s]
        j = 2 * risk
        if j > min_journey:
            continue
        j += risk_calc(s, (s[0], s[1], - s[2]))
        if j < min_journey:
            min_journey = j
    print r, len(risk_table), min_journey
    return min_journey
                        
if __name__ == '__main__':
    from time import time
    t0 = time()
    print sum(map(moon_dijkstra, [(1 << n) - 1 for n in range(1, 16)]))
    print time() - t0